GENERAL INFO ABOUT AREA AND TIMESTEP & RUN

# newarea can be 0 (no) or 1 (yes). If area is not new, it is assumed that the land outline shapefile and cropped rock outcrop is already available and doesn't need to be unzipped etc.
newarea <- 0
areaname <- "MDV"

# L8: either "Bt" or "L1"
L8downloadtype <- "Bt"


## time range parameters 
year <- c(2020:2013)
month <- c("01","02","03","04", "09", "10","11", "12")

# maximum time between satellite scenes 
timethres <- 0.6

SET PATHS

Micro IR Seagate = E IRONWOLF / ELEMENTS = D

scriptpath <- "C:/Users/mleza/OneDrive/Documents/PhD/work_packages/auto_downscaling_30m/downscale_controlscripts/data_prep/"
scriptpath_organized <- "C:/Users/mleza/OneDrive/Documents/PhD/work_packages/auto_downscaling_30m/downscale_controlscripts/data_prep/downscaling_organized/"

maindir <- "D:/downscaling_after_talk/" # this must point to an existing directory, rest is generated in setup
main <- paste0(maindir, "data_download_preprocessing/")
L8datpath <- paste0(main, "L8/")
modispath <- paste0(main, "MODIS/")
tdpath <-paste0(main, "timediff/")
cddir <- paste0(maindir, "clean_data/")
figurepath <- "C:/Users/mleza/OneDrive/Documents/PhD/work_packages/auto_downscaling_30m/paper/paper_draft/figures/"


######## set path to DEM, AOI, land outline

dempath <- "E:/new_downscaling/tiles_westcoast/" # this must point to an existing directory with dem inside
aoipath <-  "E:/new_downscaling/aoi/" # this must point to an existing directory with "Levy_MDV_actually.shp"
aoip <- list.files(aoipath, pattern="actually.shp", full.names = T)
clpath <- "E:/new_downscaling/coastline/Coastline_high_res_polygon/"

CALL SETUP

source(paste0(scriptpath_organized, "0a_setup.R"))
## Error in get(genname, envir = envir) : 
##   Objekt 'testthat_print' nicht gefunden
## OGR data source with driver: ESRI Shapefile 
## Source: "E:\new_downscaling\aoi\Levy_MDV_actually.shp", layer: "Levy_MDV_actually"
## with 1 features
## It has 1 fields
## Integer64 fields read as strings:  id 
## OGR data source with driver: ESRI Shapefile 
## Source: "E:\new_downscaling\coastline\Coastline_high_res_polygon\land_contours_MDV.shp", layer: "land_contours_MDV"
## with 14678 features
## It has 2 fields
#file.edit(paste0(scriptpath_organized, "0a_setup.R"))

Load template

# this is one raster with a complete coverage of the research area to use as a template 
template <- raster("E:/new_downscaling/clean_data/template_new.tif")

1 DEM

DEM 8m REMA

file.edit(paste0(scriptpath_organized, "1_DEM.R"))
  • Use DEM tiles “17_34_8m”, “17_35_8m”,“18_35_8m”, “18_34_8m”, “19_34_8m”, “16_35_8m”, “17_36_8m”, “19_35_8m”
  • set everything below 100m to NA
  • make 8m mosaic, everything below -50m goes
  • 3x3 mean filter
  • crop to aoi
  • fill missing values with 200m RAMP
  • resample first to 8m and cover and mask, then resample to 30m to evade sharp cuts
  • calculate slope and aspect
  • make a 20000x20000 blockmask (20km²)
dem <- raster(paste0(dempath, "DEM_30m_", areaname,"_clean_aoi_filled_mask_new.tif"))
blockmask <- raster(paste0(dempath, "blockmask_aoi.tif"))
aspect_deg <- raster(paste0(dempath, "30m_aspectMDV.tif"))
slope_deg <- raster(paste0(dempath, "30m_slopeMDV.tif"))
# aspect_rad <- raster(paste0(dempath, "30m_radians_aspectMDV.tif"))
# slope_rad <- raster(paste0(dempath, "30m_radians_slopeMDV.tif"))


mapview(aoianta, alpha.regions=0.2)+
  mapview(blockmask)+  
  mapview(aspect_deg)+
  mapview(slope_deg)+
  mapview(dem)

2 Landsat

getprocessLANDSAT() is a function(time_range)

login_USGS("MaiteLezama", "Eos300dmmmmlv")
## Login successfull. USGS ERS credentials have been saved for the current session.
login_earthdata(username="Mlezama", password="Eos300dm")
## Login successfull. NASA URS EarthData credentials have been saved for the current session.
y=1
m=2
ymid <- substring(time_range[[y]][[m]][[1]][[1]], 1, 7)
#time_range[[y]][[m]]

Landsat data are distributed in UTM - thus, the aoi is converted to UTM 57 south (WGS84). Footprints come in long lat. Per day, the footprints for “bt” data are taken from the query and their intersection with the research area is checked.

The Top of Atmosphere Brightness Temperature (toa_bt) product from Landsat 8 is used, which comes in 30m resolution, Kelvin, fill_value=“-9999” (to be rescaled by 0.1). The valid range is from 1500 to 3500.

Only those scenes are selected, where at least 30% of research area (total 20719560545m² or 20719.6km²) is covered by the scene (i.e. 6906.5km²) in the footprint. Only scenes that also show less than 20% land cloud cover are selected and written as querynew.RDS into the Landsat scene directory. Also “downloadScenes.csv” is written, which contains the date specifics on the downloaded scenes.

Time: Landsat comes in GMT and MODIS in UTC, and there is no time difference between Greenwich Mean Time and Coordinated Universal Time.

To order the scenes, a text file will be generated, with all the record ids. Those will be submitted for ordering on espa, where “Brightness Temperature - Thermal band TOA processing” and “Pixel QA” will be selected. This is a workaround, as ordering takes a lot of runtime, as a maximum of one dataset per 30 min can be ordered. Finally, all is downloaded via the package espa.tools with earthexplorer_download().

file.edit(paste0(scriptpath_organized, "2_1_selectLandsat.R"))
file.edit(paste0(scriptpath_organized, "2_2a_downloadLandsat.R"))
file.edit(paste0(scriptpath_organized, "2_2b_downloadLandsat.R"))
file.edit(paste0(scriptpath_organized, "2_3_processLandsat.R"))
file.edit(paste0(scriptpath_organized, "2_2c_place_scenes.R"))

Identify potentially good scenes:

source(paste0(scriptpath_organized, "2_1_selectLandsat.R"))

for(y in seq(year)){
  for(m in seq(month)){
    checkMyInternet()
    login_USGS("MaiteLezama", "Eos300dmmmmlv")
    try(login_earthdata(username="Mlezama", password="Eos300dm"))
    selectLANDSAT(time_range)
  }
}

“downloadScenes.csv” is where those scenes are, that passed the qualitycheck for clouds. For those scenes, I look for MODIS scenes that are located 20+/- around the L8 date. The time differences are documented in “timediff_df.csv”.

Actualized queries “MODquerymatched_msel.rds” and “L8querymatched.rds” are written out.

From the Landsat 8 query now,

Gather all record ids to submit online to espa

Download by espa ID

Only MODIS data from Aqua have a close match time.

Those 97 scenes are those from the 104 found in the second query, that had already been processed.

So, what was lost exactly in comparison to the last run, where we had > 300 scenes?

…mostly scenes from Terra (MOD) were lost, that reduces the overall time difference in the dataset. Also, the median land cloud cover decreased.

Process Landsat

Here, I’m using already downloaded images. If it would be a new download, files would be unzipped. Get date, time, filename and land cloud cover from the MODIS-matched Landsat query and write that info into timediff_df. Identify which quality assessment band codes are regarded here as clouds (cloud shadow, clouds, medium and high confidence cloud and high confidence cirrus).

Help files for selection of MODIS and documentation of time differences are written.

Brightness temperature TOA

btc_ex <- raster(paste0(L8datpath, "BTC_LC08_L1GT_057116_20190106_20190130_01_T2_bt_band10.tif"))
btc_ex_res <- raster(paste0(L8datpath, "BTC_res_LC08_L1GT_057116_20190106_20190130_01_T2_bt_band10.tif"))
LST_ex <- raster(paste0(L8datpath, "LST_LC08_L1GT_057116_20190106_20190130_01_T2_bt_band10.tif.tif"))

mapview(btc_ex, map.types="Esri.WorldImagery")+
  mapview(btc_ex_res)+
  mapview(LST_ex)
diff_BTC_LST <- LST_ex - btc_ex_res
mapview(diff_BTC_LST)
#file.edit(paste0(scriptpath_organized, "2_3_processLandsat.R"))
source(paste0(scriptpath_organized, "2_3_processLandsat.R"))

file_location <- "E:/new_downscaling/data_download_preprocessing/L8"

y=2
m=1
processLandsat(time_range, new_download=FALSE) # new_download=TRUE if all espa downloads are in one folder
workspace.size()
source(paste0(scriptpath_organized, "2a_cloudclean_Landsat_org.R"))
file.edit(paste0(scriptpath_organized, "2a_cloudclean_Landsat_org.R"))

Landsat and MODIS

MODIS LST product stems from channels 31 and 32, here is a comparison of the wavelengths covered by the thermal channels in both sensors:

MODIS

according to MODIS_LST_products_UserGuide_C5.pdf MOD11_L2, masked with the MODIS Cloud Mask data product (MOD35_L2) The LST retrieval in a MODIS swath is constrained to pixels that: (1). have nominal Level 1B radiance data in bands 31 and 32, (2). are on land or inland water, (3). are in clear-sky conditions at a confidence (defined in MOD35) of >=95% over land <= 2000m or >= 66% over land > 2000m, and at a confidence of >= 66% over lakes.

MODIS: take a look at Error_LST

use emissivity from modis product ?

compare gdal translate product with heg tool translation

CALL DOWNLOAD AND PREPROCESSING

Open questions:

#to do for final run:

\(\checkmark\) band 10/11 reference (literature)